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Abstract 

Starting from an atomistic approach we have de- 
rived a hierarchy of successively more simplified con- 
tinuum elasticity descriptions for modeling the me- 
chanical properties of suspended graphene sheets. 
The descriptions are validated by applying them 
to square graphene-based resonators with clamped 
edges and studying numerically their mechanical re- 
sponses. Both static and dynamic responses are 
treated. We find that already for deflections of the 
order of 0.5A a theory that correctly accounts for 
nonlinearities is necessary and that for many pur- 
poses a set of coupled Duffing-type equations may be 
used to accurately describe the dynamics of graphene 
membranes. 

Introduction. Recent progresses in fabricating 
graphene, a monolayer of graphite, have considerably 
boosted the attention for this material [T| . Among its 
unique features are remarkable electronic properties 
[3J [3] , which make graphene of considerable interest 
for both fundamental science and technological ap- 
plications. Moreover, its exceptionally large mechan- 
ical strength [4j and ability to sustain large electrical 
currents can be of great value in the broad field of 
nanodeviccs. In particular, in the field of nanoelec- 
tromechanical systems (NEMS), graphene-based me- 
chanical resonators were recently demonstrated [H [5] 
and theoretical work indicates a strong coupling be- 
tween deformation and intrinsic electronic properties 
of graphene [7] . 

Most research on graphene has hitherto focused on 



the electronic properties of graphene, and less at- 
tention has been directed to mechanical properties. 
For modeling NEMS a reliable and efficient descrip- 
tion of the mechanical response of nanocarbons to 
external forces is essential [H [9]. While continuum 
elasticity theory has been applied succesfully to the 
study of mechanical properties of nanotubes for a 
long time, [lOl [HI H2] it has only recently been ap- 
plied to graphene membranes [T31 [TH Q21 [TB] . 

In this work we first formulate a general nonlinear 
elasticity theory for graphene sheets starting from an 
atomistic description in terms of a valence force field 
model. Through successive approximations, we then 
derive simplified models which are easier to solve and 
to study analytically. We then apply these continuum 
descriptions to a drum-like resonator and investigate 
the mechanical response in both the static and dy- 
namic cases. The results obtained agrees well with 
the experimentally observed responses [TTJ [T5] . 
Continuum elasticity model for graphene. 
The interaction potential between carbon atoms in 
graphene can be modeled by a valence force field 
model [19l [20] where the potential energy U sp 2 be- 
tween sp 2 bonds is given by 
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Here the index i labels the carbon atoms while in- 
dices j and k are bond labels for the nearest neigh- 
bor atoms of i. Thus, fij, j — 1,2,3, are the three 
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bond vectors that connect the atom i to its nearest 
neighbors. The parameters a, and 7 are constants, 
N at the number of atoms, ao = 1.421 A is the equi- 
librium bond length in graphite, and D L — ._ 1 fy 
is the dangling bond vector. The first two terms of 
Eq. (fTJ) represent the energy cost necessary to change 
the length and angle between covalent C-C bonds. 
In the continuum description, we refer to it as the 
stretching energy. The last term of Eq. ([T|) gives the 
energy cost necessary to change the angle between 
p 2 -orbitals, which are approximately normal to the 
graphene surface. We refer to energy contributions 
from this term as bending energy. 

Provided that the length scale of the deforma- 
tion is large compared to the lattice spacing (long 
wavelength limit), continuum theory can be used for 
graphene [15] and we can write 

U sp 2 = J dxdy W [u(x,y)] 

by parametrizing the deformed surface as u(x,y) = 
[u(x, y),v(x, y),w(x,y)]. The elastic energy density 
Wo can be divided into stretching and bending con- 
tributions, Wo = Wq + W(f . The stretching energy 
density can be written as [21j 

Eh 

W ° S = 2(1 - v*) i E L+E 2 yy +2^E xx E yy +2(l-iy)El y ], 

(2) 

where the components of the Green strain tensor Eij 
are 

E xx = u x + {u 2 x +v 2 x +w 2 x )/2, 

E xy = (u y + v x + u x u y + v x v y + w x w y )/2, 

Eyy = V y + (u 2 y +V 2 y +W 2 y )/2. 

Here the subscripts on u, v, w denote differentiation, 
i.e. u x = du/dx etc., and the coefficient Eh rep- 
resents the Young modulus multiplied by the thick- 
ness of thin plate theory [22l [23] . It is important 
to note that in our theory Eh is a single parameter 
and we do not consider any thickness in our formu- 
lation. The Poisson ratio is denoted by v. Both Eh 
and v are related to the Lame parameters A and p 
as Eh = 4p(X + /i)/(A + 2p) and v = A/(A + 2p). 
The Lame parameters are determined from a and (3 



as 6p = V3(a + 8/3) and 6A = V3(a - 4/3) (for de- 
tails, see [21]). Using values from Ref.[T9] (a = 155.9 
J/m 2 and (3 = 25.5 J/m 2 ), gives fx = 103.89 J/m 2 
and A = 15.55 J/m 2 . The continuum theory is in an 
excellent agreement with molecular dynamics simu- 
lations of graphene [15] . 

The bending energy density can be approximated 

as 

W B = ^2a H) 2 , 

where k — ^/2>aff)/2 w 0.8 eV is the bending ridigity 
and H = (w xx + w yy )/2 is the local mean curvature. 

In clamped graphene-based NEMS applications 
(e.g: mechanical resonators), we find that the stretch- 
ing energy is much greater than the bending energy. 
In thin plate theory, the bending regime (linear the- 
ory) is valid only for out-of-plane deflections less than 
the plate thickness [23]. Following this argument, for 
graphene the linear regime is almost nonexistent and 
the nonlinear stretching regime is dominant. We have 
determined this to be the case for deflections in ex- 
cess of 0.5 A. Thus, in the following we neglect the 
bending contribution. 

For a given applied external body force Fq, the 
equilibrium shape can be obtained by direct mini- 
mization of the free energy functional 

!F[u(x, y)] = J dxdyWo — J dxdypQm^ 1 F ■ u, 

where po is the mass density and m c is the carbon 
mass. Alternatively, the equilibrium shape may be 
found by solving for the dynamics of the system in- 
cluding dissipation. We will here work with the lat- 
ter approach. The dynamic equation in the stretch- 
ing regime, where Eq. @ is used as the elastic en- 
ergy density, can be obtained from standard varia- 
tional principles. The corresponding equation of mo- 
tion, with an added phcnomenological damping term 
cu(x,y), is 

u(x,y) + cu(x,y) = pQ lr DP[u(x,y)}+m c ~ 1 F (x,y,t), 

(3) 

where P is the Piola stress tensor. In cartesian coor- 
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dinates (xy), the Piola stress tensor is 

PXX = (I + U X )[(\ + 2[l)E XX + XEyy] +2(lUyE X y, 

P X y = 2/i(l + Vy)E Xy + V X [(X + 2fJ,)E XX + XEyy], 

PxZ = XW X (E XX + Eyy) + 2[IW X E XX + 2/IWyExy, 
Py X = 2fJ,(l + U X )E X y+Uy[(X + 2fJ,)Eyy + \E xx ], 

Pyy = (l + V y )[(X + 2(J,)E y y + XE xx ]+2(J,V x E x y, 

Py Z = XWy(E xx + Eyy) +2fJ,WyEyy +2flW X E x 0) 

and the linear differential operator V acts on P as 



VP = 



(d x P XX +dyPy x )x 



Dirichlet boundary conditions can be imposed 
through u = uo(t) and Neuman boundary conditions 
through no ■ P = Po(t), where uo(t) and Po(t) are 
specified on the boundary. 

We refer to Eqs. ([3]) and (fj| as the full non- 
linear system or general elasticity equations. Since 
the general expresions for the Piola stress tensor are 
quite cumbersome and make the equations difficult to 
solve. The full equations can be approximated by a 
simpler set of equations where we neglect the second 
order contributions of the in-plane displacements u 
and v. This yields the following expressions for the 
components of the strain tensor 
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u x + w x /2, 

(Uy +V X + W x Wy)/2, 
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and the Piola tensor 
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(A + 2[l)E XX + XEyy, 
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P XZ ~ XW X {E XX + Eyy) + 2fJb{W X E XX + WyE X y), 

Pyx ~ ^f^E X y , 



p 



^E xx + (A + 2fJ,)Eyy, 



XWy{E : 



1/ 1 . T Eyy ) 



2fl(WyEyy + W X E Xy \6) 



The resulting equations of motion for u, v and w are 
known in thin plate theory [22J, [23] as the nonlinear 
von Karman equations. It is worth mentioning that 



we have arrived at these equations without treating 
the graphene as a thin plate with some thickness. 

We may further simplify the above expressions by 
completely removing the in-plane displacements. The 
resulting nonlinear equation for w(x,y,t) is given by 



)(x,y,t) + cw(x,y,t) - p 1 d x (w x T x ) 



(7) 



where 

T x = {X + 2 f i)5 x + X5 y + (X/2 + [i)(w 2 x +w y ), 
T y = (X + 2 f j,)Sy + XS x + (X/2 + n)(wl+wl). 

Here we have introduced the constants 6 X and S y rep- 
resenting initial strains in the x and y directions, re- 
spectively. Such strains may appear during the man- 
ufacturing process of graphene resonators. The func- 
tions T x and T y are tensions in the x and y directions 
induced by stretching of the graphene. We refer to 
Eq. as the out-of-plane approximation. 

If one mode is expected to dominate the out-of- 
plane deformations, the out-of-plane approximation 
may be projected onto this mode to obtain an ordi- 
nary rather than a partial differential equation. If 
the dominant mode is the fundamental, we can write 
for a square sheet (—a < x, y, < a) 

w{x, y, t) = w{t) cos( — ) cos — . 

2a 2a 

By using this Ansatz, we obtain a Duffing equation 
for the amplitude of the fundamental mode, 



w(i) + cw(t) + uj Q w(t) 



where 



5tt 4 (A + 2^) 3 F{t) 



128a 4 po 



w (*) 
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F W = \ II F(x,y,t)cos(^)cos{^)dxdy 
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is the overlap of the driving force with the shape of 
the fundamental mode, and the resonant frequency 
luq is given by 
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7T (\ + fi)(6 + 5 2 /2) 
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and we have assumed S x = 5 y = 6. Notice if sev- 
eral modes are included in the Ansatz, we obtain a 
set of coupled Duffing-type equations for the mode 
amplitudes. 

Numerical results: Statics. We now proceed to 
investigate the accuracies of the different approx- 
imations — full nonlinear system ([3]|4]), von Kar- 
man approximation (JSHHJ), out-of-plane approxima- 
tion |(7J), and the Duffing equation (jHJ). For illustra- 
tion, we consider in the simulations a clamped square 
graphene sheet of side 2a = 1 /im, (Dirichlet bound- 
ary condition with Uo(t) = on all edges) subject to 
uniform pressure F dc z/A where A — 3\/3ao/4 is the 
area per carbon atom in graphene. In all computa- 
tions we consider a pre-tensile graphene with initial 
strains of 5 X = S y = 0.5%. We solve the Eq. ^ nu- 
merically by using the spectral method 24. 25j with 
81 global basis functions. 

The tension in the graphene sheet varies depend- 
ing on how much stretching is imposed by the exter- 
nal forces. For pre-tensile graphene, small deflections 
do not change the built-in tension and the graphcnc- 
based resonator behaves as a linear tensile elastic 
membrane [26 . However, when the deflections are 
large enough, the stress changes significantly from 
the initial value. To illustrate this, we present Figs. 
[T] and [2] The first figure shows the stress distribu- 
tions for two values of driving forces F dc per carbon 
atom, Fdc = 1 fN and Fd c = 1 pN. For the smaller 
force, the deflection at equilibrium is small and the 
tension P xx is almost uniform and equal to the ini- 
tial value 1.2 N/m. This agrees with the theoretical 
value of a = Eh/{l-v)5 w 1.2 N/m, where 8 = 0.5% 
is the initial strain of the pre-tensile graphene. On 
the other hand, for the larger force Fd c = 1 pN, the 
equilibrium built-in tension P xx , as well as the other 
components of the Piola stress tensor, is not uniform 
and it varies from 2 N/m to 14 N/m. The inset in Fig. 
[2] shows the dependence between the average value of 
the tension P xx and the driving force Fd c . Note that 
for small F<j c the tension is almost constant and equal 
to 1.2 N/m, but for large Fd c the tension varies as Fd c 
as P xx ~ fJ c . 

Next, we present the results for the vertical deflec- 
tion at the center of the graphene sheet as function of 
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Figure 1: Tension distributions for small deflections 
with Fdc = 1 fN per carbon atom (upper plots) and 
large deflections with Fdc = 1 pN (lower plots). Note 
the very different colour scales. 

the external force Fd c - We find that for a pre-tensile 
graphene sheet the response is linear for small values 
of Fdc- This is due to an almost uniform tension, in- 
dependent of the amount of deflection, see Fig. [2] On 
the other hand, when the out-of-plane deflections are 

1/3 

large, we see a nonlinear response w ~ F,' , where 
w is the deflection at the center of the graphene. 

Finally, evaluations of the bending energy for the 
equilibrium shapes, for a range of Fdc between 0.1 fN 
and 1 pN, show that the stretching energy is at least 
four orders of magnitude larger than the bending en- 
ergy. This validates the assumption that bending en- 
ergy in clamped graphene-based nanoresonators may 
be neglected. In addition, we find that all the sim- 
plified mechanical descriptions show good agreement 
with the general theory in the static case. 
Numerical results: Dynamics. In this section, we 
study the dynamic response of pre-tensile clamped 
graphene subject to small and large deformations. 
For small deformations, when the tension is set by the 
initial strain, the clamped graphene sheet behaves as 
a linear tensile elastic membrane [26] . On the other 
hand, when the out-of-plane deformations are large 
and tension depends on the amount of deflection, a 
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Figure 2: Elastostatics study for clamped pre-tensile 
graphene. The figure shows the vertical deflection 
at the center of the graphene sheet in response to 
a uniform force ranging from 0.1 fN to 1 pN per 
atom. Good agreement is obtained between the sim- 
plified models and the full nonlinear system. The 
inset shows the dependence between the tension on 
the driving force. 

nonlinear response governs the dynamic behaviour. 
Here we show that the dynamics exhibits Duffing- 
type behavior, e.g. we observe the characteristic in- 
stability in the amplitude when we sweep the driving 
frequency. 

We divide the analysis into three parts. First, we 
consider a driving force with the same spatial distri- 
bution as the fundamental mode. Here we have only 
one mode present which makes this case easy to anal- 
yse. In the second part we consider a uniform force 
where several modes couple to the external force. 
Here we observe the presence of up to three modes 
with similar Dufnng-type response as the fundamen- 
tal mode. Finally, we study small oscillations around 
an equilibrium determined by the dc-component of 
the exciting forces. We evaluate the resonance fre- 
quencies as function of the Fd c and the amplitude of 
oscillation as a function of the ac-component of the 
driving force F ac . 

For a driving force with the same spatial depen- 
dence as the fundamental mode, 

F(x,y,t)=F(t)cos(^)cos(^), 
la la 
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Figure 3: Amplitude of the fundamental mode for 
different driving forces strength. Observe the ampli- 
tude instability for large values of F ac . All points 
were obtained starting from a configuration with the 
sheet at rest. If the dynamic state of the system is 
kept between frequency steps a hysteretic response 
obtains. 

the out-of-plane and Duffing approximations are 
equivalent. We perform numerical computations for 
small and large harmonic driving forces (F(t) = 
F ac cos(wi)) . The results are shown in Fig. [3l The 
figure shows that for small driving force F ac — 0.5 fN, 
the response is that of a linear harmonic resonator. 
However, for the larger driving forces of F ac = 6 fN 
and 20 fN, we see a nonlinear response which mani- 
fests as an amplitude instability at a critical, ampli- 
tude dependent frequency; e.g., at / = 1100 MHz for 
F ac = 6 fN. The Duffing equation predicts this insta- 
bility at a slightly higher frequency. This is because 
the tension is slightly reduced when the in-plane dis- 
placements are included (general elasticity theory or 
von Karman approximation), see Fig. O There is 
an excellent agreement between the general elasticity 
theory and the von Karman approximation. 

Next, we study the case of a uniform driving force. 
Here several modes couple to the driving forces. The 
fundamental mode has the greatest overlap followed 
by the modes 

.3irx . .ny ,irx . ,37ra. 
cos (^— )cos( — ),cos( — )cos( — ), 
la la la la 
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and 

3irx . , 3ny 
C0S( ^ )COS( ^ ) - 

These modes are found by plotting the stretching en- 
ergy against the driving frequency. This is shown in 
Fig. [U The figure shows three peaks correspond- 
ing to the different modes: the peak at 1200 MHz 
belongs to the fundamental mode, and the peaks at 
2360 MHz and 2840 MHz correspond to the other 
modes mentioned above. Notice that the first peak 
at 350 MHz does not correspond to any mode but to 
a i sub-harmonic of the fundamental. The presence 
of this sub-harmonic is due to the cubic nonlinear 
term in Eq. ©. 

Finally, we study small oscillations around an equi- 
librium shape determined by a dc-part of the driving 
force Fd c + F ac cos(wt), typical for a NEMS resonator. 
We assume a spatially uniform driving force. We find 
that the resonance frequency for small oscillations in- 
creases as we increase the stretching or deflection at 
the equilibrium position. This may be used to tune 
the frequency of the oscillator by changing the dc-bias 

f pl/3 
aS Jres OC " dc ' 

Next, we study the dependence between the am- 
plitude of the oscillations as function of the driving 
ac- force F ac . We find that there is a linear response 

1/3 

for small F ac and then it tends to follow the law F a i . 
This is shown in Fig. [5] Here the driving frequency 
is the resonance frequency for Fd c — 250 fN, which is 



3360MHZ. 

Conclusions. We have derived and analyzed a 
nonlinear finite elasticity theory for graphene res- 
onators, both for elastostatics and elastodynamics 
problems. Moreover, we have studied how this gen- 
eral elasticity theory can be simplified to more easily 
solvable equations. In particular, the out-of-plane ap- 
proximation, Eq. ([7]), gives good agreement with the 
general elasticity theory while maintaining the ad- 
vantage of being computationally efficient. We have 
also shown that the dynamic response of clamped 
graphene resembles that of coupled Duffing-type res- 
onators. 
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